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<H Summary. In this paper we test a special-relativistic formulation of Smoothed 

j"*^ Particle Hydrodynamics (SPH) that has been derived from the Lagrangian of an 

^__( ideal fluid. Apart from its symmetry in the particle indices, the new formulation 

differs from earlier approaches in its artificial viscosity and in the use of special- 

I , ' relativistic "grad-h-terms" . In this paper we benchmark the scheme in a number of 

^>j demanding test problems. Maybe not too surprising for such a Lagrangian scheme, 

KH it performs close to perfectly in pure advection tests. What is more, the method 

(-H produces accurate results even in highly relativistic shock problems. 
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1 Introduction 

> 

Qs^ Relativity is a crucial ingredient in a variety of astrophysical phenomena. For 

f^ example the jets that are expelled from the cores of active galaxies reach 

^^O velocities tantalizingly close to the speed of light, and motion near a black 

" . hole is heavily influenced by space-time curvature effects. In the recent past, 

VO substantial progress has been made in the development of numerical tools 

^^ to tackle relativistic gas dynamics problems, both on the special- and the 

general-relativistic side, for reviews see [IHllIlllS]- Most work on numerical 
relativistic gas dynamics has been performed in an Eulerian framework, a 
couple of Lagrangian smooth particle hydrodynamics (SPH) approaches do 
exist though. 
^ In astrophysics, the SPH method has been very successful, mainly because 

^ of its excellent conservation properties, its natural flexibility and robustness. 

Moreover, its physically intuitive formulation has enabled the inclusion of 
various physical processes beyond gas dynamics so that many challenging 
multi-physics problems could be tackled. For recent reviews of the method 
we refer to the literature |2ll El] ■ Relativistic versions of the SPH method 
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were first applied to special relativity and to gas flows evolving in a fixed 
background metric [M HSl [H El HIST] . More recently, SPH has also been 
used in combination with approximative schemes to dynamically evolve space- 
time [fl [TOl [HI m [23 H H I3j • 

In this paper we briefly summarize the main equations of a new, special- 
relativistic SPH formulation that has been derived from the Lagrangian of an 
ideal fluid. Since the details of the derivation have been outlined elsewhere, 
we focus here on a set of numerical benchmark tests that complement those 
shown in the original paper [5B]. Some of them are "standard" and often 
used to demonstrate or compare code performance, but most of them are 
more violent — and therefore more challenging — versions of widespread test 
problems. 



2 Relativistic SPH equations from a variational principle 

An elegant approach to derive relativistic SPH equations based on the dis- 
cretized Lagrangian of a perfect fluid was suggested in [25] . We have recently 
extended this approach 128, 29^ by including the relativistic generalizations of 
what are called "grad-h-terms" in non-relativistic SPH [351 [53] ■ For details of 
the derivation we refer to the original paper [28j and a recent review on the 
Smooth Particle Hydrodynamics method [271 . 

In the following, we assume a flat space-time metric with signature (-,+,+,+) 
and use units in which the speed of light is equal to unity, c = 1. We re- 
serve Greek letters for space-time indices from 0...3 with being the temporal 
component, while i and j refer to spatial components and SPH particles are 
labeled by a, b and k. 

Using the Einstein sum convention the Lagrangian of a special-relativistic 
perfect fluid can be written as [13] 

Lpt^sr = - f T^'U^.U, dV, (1) 

where 

T^"" = (n[l + u{n, s)] + Pp^'U" + Prf (2) 

denotes the energy momentum tensor, n is the baryon number density, u 
is the thermal energy per baryon, s the specific entropy, P the pressure and 
If^ — dx^ /dr is the four velocity with r being proper time. All fluid quantities 
are measured in the local rest frame, energies are measured in units of the 
baryon rest mass energjFl mgc^. For practical simulations we give up general 
covariance and perform the calculations in a chosen "computing frame" (CF). 
In the general case, a fluid element moves with respect to this frame, therefore. 



^ The appropriate mass mo obviously depends on the ratio of neutrons to protons, 
i.e. on the nuclear composition of the considered fluid. 
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the baryon number density in the CF, N, is related to the local fluid rest frame 
via a Lorentz contraction 

N = 7n, (3) 

where 7 is the Lorentz factor of the fluid element as measured in the CF. The 
simulation volume in the CF can be subdivided into volume elements such that 
each element b contains 1^1, baryons and these volume elements, AVb = h't/Nf,, 
can be used in the SPH discretization process of a quantity /: 

b ^ 

where the index labels quantities at the position of particle 5, r;,. Our notation 
does not distinguish between the approximated values (the / on the LHS) and 
the values at the particle positions (/f, on the RHS). The quantity h is the 
smoothing length that characterizes the width of the smoothing kernel W , for 
which we apply the cubic spline kernel that is commonly used in SPH [531 1211 • 
Applied to the baryon number density in the CF at the position of particle a, 
Eq. Q yields: 

Na^N{ra) = Y,^hW{\ra~rklK). (5) 

6 

This equation takes over the role of the usual density summation of non- 
relativistic SPH, p{ra) = X]fc"^fc^(ka ~ rb\,h). Since we keep the baryon 
numbers associated with each SPH particle, i^b, fixed, there is no need to 
evolve a continuity equation and baryon number is conserved by construction. 
If desired, the continuity equation can be solved though, see e.g. [3]. Note that 
we have used a's own smoothing length in evaluating the kernel in Eq. ([5]). 
To fully exploit the natural adaptivity of a particle method, we adapt the 
smoothing length according to 

where 77 is a suitably chosen numerical constant, usually in the range between 
1.3 and 1.5, and D is the number of spatial dimensions. Hence, similar to the 
non-relativistic case [32 [SJ , the density and the smoothing length mutually 
depend on each other and a self-consistent solution for both can be obtained 
by performing an iteration until convergence is reached. 
With these prerequisites at hand, the fluid Lagrangian can be discretized 

[23 [27] 

isPH.sr = - V" — [1 +u(n6,Sb)]. (7) 

Using the first law of thermodynamics one finds (for a detailed derivation see 
Sec. 4 in [271) for the canonical momentum per baryon 

Sa= ?; — = laVa 1 + Ua H , (8) 
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which is the quantity that we evolve numerically. Its evolution equation follows 
from the Euler-Lagrange equations, 

d dL dL 
dt dva dra 

as [27] 

f^ - E - (]^v.iy..(..) + ^^^.wm) , (10) 

where the "grad-h" correction factor 

dhb ^^ dWbkjhb) 
dNb ^ dhb 

was introduced. As numerical energy variable we use the canonical energy per 
baryon, 

Ca = 7a 1 + Ua + ] - ^:r = '" ^^ ' ^ a A (12) 

V naj Na 7a 

which evolves according to [27] 

1^ - E - (^ ■ v.«'..('..) + ^ ■ -^Mk,)) . (13) 

As in grid-based approaches, at each time step a conversion between the nu- 
merical and the physical variables is required [H [55] . 

The set of equations needs to be closed by an equation of state. In all of the 
following tests, we use a polytropic equation of state, P = {r—l)nu, where F 
is the polytropic exponent (keep in mind our convention of measuring energies 
in units of toqc^). 



3 Artificial dissipation 

To handle shocks, additional artificial dissipation terms need to be included. 
We use terms similar to [4] 

-^] ^-y^^bHabyaWab with 77,fc = - -=^ (S: - 5^ ) • 6,6 (14) 

dt J diss b ^''b 

and 

^) ^-E^^*'""-^^^^ ^itl^ '^ab^-^iC-eDiab. (15) 

"'t J diss b "■'' 
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Here K is a numerical constant of order unity, Wsjg an appropriately chosen 
signal velocity, see below, Nab = (^a + ^^6)/2, and iab = (^a - f f,)/ka - f b| 
is the unit vector pointing from particle b to particle a. For the symmetrized 
kernel gradient we use 

1 



"^aWab = 2 l^aWabiha) + "^ aWabih)] ■ (16) 



Note that in [4J V aWab{hab) was used instead of our VaWab, in practice we 
find the differences between the two symmetrizations negligible. The stars at 



the variables in Eqs. (14 1 and (151 indicate that the projected Lorentz factors 



it - /^ / . ,, (17) 

^yl - [Vk ■ Bab)'' 

are used instead of the normal Lorentz factor. This projection onto the line 
connecting particle a and b has been chosen to guarantee that the viscous 
dissipation is positive definite [1]. 

The signal velocity, Wgig, is an estimate for the speed of approach of a signal 
sent from particle a to particle 6. The idea is to have a robust estimate that 
does not require much computational effort. We use [55] 

Vsig,ab = max(aa,a6), (18) 

where 

.± ^n _L\± 



a. 



max(0,±A^) (19) 



with Aj. being the extreme local eigenvalues of the Euler equations 



^± __ Vk ± Cs,fc 
1 ± VkCsh 



A± = ^', l'" (20) 



and Cs,fc being the relativistic sound velocity of particle k. These ID estimates 
can be generalized to higher spatial dimensions, see e.g. [201. The results 
are not particularly sensitive to the exact form of the signal velocity, but in 



experiments we find that Eq. ( 18 ) yields somewhat crisper shock fronts and 
less smeared contact discontinuities (for the same value of K) than earlier 
suggestions [1]. 

Since we are aiming at solving the relativistic evolution equations of an ideal 
fluid, we want dissipation only where it is really needed, i.e. near shocks where 
entropy needs to be producecPl To this end, we assign an individual value of 
the parameter K to each SPH particle and integrate an additional differential 
equation to determine its value. For the details of the time-dependent viscosity 
parameter treatment we refer to |28j . 



^ A description of the general reasoning behind artificial viscosity can be found, for 
example, in Sec. 2.7 of [27] 
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4 Test bench 

In the following we demonstrate the performance of the above described 
scheme at a slew of benchmark tests. The exact solutions of the Riemann 
problems have been obtained by help of the RIEMANN_VT.f code provided 
by Marti and Miiller f20|. Unless mentioned otherwise, approximately 3000 
particles are shown. 

4.1 Test 1: Riemann problem 1 

This moderately relativistic (maximum Lorcntz factor 7niax ~ 1-4) shock 
tube has become a standard touch-stone for relativistic hydrodynamics codes 
[13 mi H EOl 13 I2D]- It uses a polytropic equation of state (EOS) with an 
exponent of -T = 5/3 and [P, N,v]i, — [40/3, 10, 0] for the left-hand state and 
[P, N, u]pj = [10~^, 1, 0] for the right-hand state. As shown in Fig. [I] the nu- 
merical solution at t = 0.35 (circles) agrees nearly perfectly with the exact 
one. Note in particular the absence of any spikes in u and P at the contact 
discontinuity (near x ~ 0.25), such spikes had plagued many earlier relativis- 
tic SPH formulations [17l [31] . The only places where we see possibly room for 
improvement is the contact discontinuity which is slightly smeared out and 
the slight over-/undershoots at the edges of the rarefaction fan. 
In order to monitor how the error in the numerical solution decreases as a 
function of increased resolution, we calculate 

Li^^^Y.\''b-vUrh)l (21) 



part 



b 



where iVpart is the number of SPH-particles, Vh the (ID) velocity of SPH- 
particle h and Ve-^{ri,) the exact solution for the velocity at position r;,. The 
results for Li are displayed in Fig. 2 The error Li decreases close to ex A^" ^ 
(actually, the best fit is Li ex Apg'J^^, which is what is also found for Eulerian 
methods in tests that involve shocks. Therefore, for problems that involve 
shocks we consider the method first-order accurate. The order of the method 
for smooth flows will be determined in the context of test 6. 

4.2 Test 2: Riemann problem 2 

This test is a more violent version of test 1 in which we increase the initial left 
side pressure by a factor of 100, but leave the other properties, in particular 
the right-hand state, unchanged: [F, p, uJl = [4000/3,10,0] and [P,p,v]^ = 
[10~^, 1,0]. This represents a challenging test since the post-shock density is 
compressed into a very narrow "spike" , at i = 0.35 near x sa 0.35. A maximum 
Lorentz-factor of 7niax ~ 3.85 is reached in this test. 
In Fig. [3] we show the SPH results (circles) of velocity u, specific energy u, the 
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Figure 1. Results of the relativistic shock tube of test 1 at i = 0.35: SPH results 
(circles) vs. exact solution (red line). From left to right, top to bottom: velocity (in 
units of c), specific energy, computing frame baryon number density and pressure. 



computing frame number density N and the pressure P at t = 0.35 together 
with the exact solution of the problem (red line) . Again the numerical solution 
is in excellent agreement with the exact one, only in the specific energy near 
the contact discontinuity occurs some smearing. 



4.3 Test 3: Riemann problem 3 

This test is an even more violent version of the previous tests. We now increase 
the initial left side pressure by a factor of 1000 with respect to test 1, but leave 
the other properties unchanged: [P, p, u]l = [40000/3,10,0] and [P,p,v]b. = 
[10~^,1,0]. The post-shock density is now compressed into a very narrow 
"needle" with a width of only « 0.002, the maximum Lorentz factor is 6.65. 
Fig. H shows the SPH results (circles) of velocity v, specific energy u, the 
computing frame number density N and the pressure P at t — 0.2 together 
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Figure 2. Decrease of the error as defined in Eq. (21 1 as a function of particle 



number for the relativistic shock tested in Riemann problem 1. The error decreases 
close to Li ex iVp-Jrt. 



with the exact solution (red hne). The overaU performance in this extremely 
challenging test is still very good. The peak velocity plateau with v « 0.99 
(panel 1) is very well captured, practically no oscillations behind the shock 
are visible. Of course, the "needle-like" appearance of the compressed density 
shell (panel 3) poses a serious problem to every numerical scheme at finite 
resolution. At the applied resolution, the numerical peak value of N is only 
about half of the exact solution. Moreover, this extremely demanding test 
reveals an artifact of our scheme: the shock front is propagating at slightly 
too large a speed. This problem decreases with increasing numerical resolution 



and experimenting with the parameter K of Eqs. ( 14 1 and ( 15 1 shows that it is 
related to the form of artificial viscosity, smaller offsets occur for lower values 
of the viscosity parameter K. Here further improvements would be desirable. 



4.4 Test 4: Sinusoidally perturbed Riemann problem 

This is a more extreme version of the test suggested by |^ . It starts from an 
initial setup similar to a normal Riemann problem, but with the right state 
being sinusoidally perturbed. What makes this test challenging is that the 
smooth structure (sine wave) needs to be transported across the shock, i.e. 
kinetic energy needs to be dissipated into heat to avoid spurious post-shock 
oscillations, but not too much since otherwise the (physical!) sine oscillations 
in the post-shock state are not accurately captured. We use a polytropic ex- 
ponent of _r = 5/3 and 

[P, TV, w]L = [1000,5,0] and [P,iV, w]^ = [5, 2 + 0.3 sin(502:), 0]. (22) 

as initial conditions, i.e. we have increased the initial left pressure by a factor 
of 200 in comparison to f6l. The numerical result (circles) is shown in Fig. [5] 
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Figure 3. Same as previous test, but tlie initial left hand side pressure has been 
increased by a factor of 100. SPH results (at t = 0.35) are shown as circles, the 
exact solution as red line. From left to right, top to bottom: velocity (in units of c), 
specific energy, computing frame baryon number density and pressure. 

together with two exact solutions, for the right-hand side densities N^ = 2.3 
(solid blue) and N-^ = 1.7 (solid red). All the transitions are located at the 
correct positions, in the post-shock density shell the solution nicely oscillates 
between the extremes indicated by the solid lines. 



4.5 Test 5: Relativistic Einfeldt rarefaction test 

The initial conditions of the Einfeldt rarefaction test [7j do not exhibit dis- 
continuities in density or pressure, but the two halfs of the computational 
domain move in opposite directions and thereby create a very low-density re- 
gion around the initial velocity discontinuity. This low-density region poses a 
serious challenge for some iterative Riemann solvers, which can return negative 
density/pressure values in this region. Here we generalize the test to a relativis- 
tic problem in which left/right states move with velocity -0.9/-f 0.9 away from 
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Figure 4. Same as first shock tube test, but the initial left hand side pressure has 
been increased by a factor of 1000. SPH results (at t — 0.2) are shown as circles, the 
exact solution as red line. From left to right, top to bottom: velocity (in units of c), 
specific energy, computing frame baryon number density and pressure. 



the central position. For the left and right state we use [P, n, w]l = [1, 1, —0.9] 
and [P, n, fjfj, = [1, 1,0.9] and an adiabatic exponent oi F = 4/3. Note that 
here we have specified the local rest frame density, n, which is related to the 
computing frame density by Eq. (Is]). The SPH solution at i = 0.2 is shown in 
Fig. [6] as circles, the exact solution is indicated by the solid red line. Small os- 
cillations are visible near the center, mainly in v and u, and over-/undershoots 
occur near the edges of the rarefaction fan, but overall the numerical solution 
is very close to the analytical one. In its current form, the code can stably 
handle velocities up to 0.99999, i.e. Lorentz factors 7 > 200, but at late times 
there are practically no more particles in the center (SPH's approximation to 
the emerging near- vacuum) , so that it becomes increasingly difficult to resolve 
the central velocity plateau. 
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Figure 5. Riemann problem where the right-hand side is periodically perturbed. 
The SPH solution is shown as circles, the exact solution for Riemann problems with 
constant RHS densities A''r — 2.3 (blue) and A^r = 1.7 (red) are overlaid as solid 
lines. 



4.6 Test 6: Ultra-relativistic advection 



In this test problem we explore the ability to accurately advect a smooth den- 
sity pattern at an ultra-relativistic velocity across a periodic box. Since this 
test does not involve shocks we do not apply any artificial dissipation. We use 
only 500 equidistantly placed particles in the interval [0, 1], enforce periodic 
boundary conditions and use a polytropic exponent of T = 4/3. We impose 
a computing frame number density N{x) — Nq + ^ sin(27rx) -I- ^ sin(47ra:), a 
constant velocity as large as u = 0.99999999, corresponding to a Lorentz 
factor of 7 sa 7071, and instantiate a constant pressure corresponding to 
Pq = (r — l)noUo, where tiq = Nq/j and A'o = 1 and uq = 1. The spe- 
cific energies are chosen so that each particle has the same pressure Pq. With 
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Figure 6. Relativistic version of tlie Einfeldt rarefaction test. Initially the flow has 
constant values of n = 1, P = 1 everywhere, vl ~ —0.9 and vr = 0.9. 



these initial conditions the specified density pattern should just be advected 
across the box without being changed in shape. 

The numerical result after 50 times (blue circles) and 100 times (green tri- 
angles) crossing the interval is displayed in Fig. [7J left panel. The advection 
is essentially perfect, no deviation from the initial condition (solid, red line) 
is visible. 

We use this test to measure the convergence of the method in the case of 
smooth flow (for the case involving shocks, see the discussion at the end of 
test 1). Since for this test the velocity is constant everywhere, we use the com- 
puting frame number density N to calculate Li similar to Eq. (21). We find 
that the error decreases very close to Li (x iV~^, see Fig.JTJ right panel, which 
is the behavior that is theoretically expected for smooth functions, the used 
kernel and perfectly distributed particles [12] (actually, we find as a best-fit 
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Figure 7. Left: Ultra-relativistic advection (v = 0.99999999, Lorentz factor 7 = 
7071) of a density pattern across a periodic box. The advection is essentially perfect, 
the patterns after 50 (blue circles) and 100 (green triangles) times crossing the box 
are virtually identical to the initial condition (red line). Right: Decrease of the 
Li error as a function of resolution, for smooth flows the method is second-order 
accurate. 



exponent -2.07). Therefore, we consider the method second-order accurate for 
smooth flows. 



5 Conclusions 



We have summarized a new special-relativistic SPH formulation that is de- 
rived from the Lagrangian of an ideal fluid |28j . As numerical variables it uses 
the canonical energy and momentum per baryon whose evolution equations 
follow stringently from the Euler-Lagrange equations. We have further applied 
the special-relativistic generalizations of the so-called "grad-h-terms" and a 
refined artificial viscosity scheme with time dependent parameters. 
The main focus of this paper is the presentation of a set of challenging bench- 
mark tests that complement those of the original paper [28 . They show the 
excellent advection properties of the method, but also its ability to accurately 
handle even very strong relativistic shocks. In the extreme shock tube test 3, 
where the post-shock density shell is compressed into a width of only 0.1 % 
of the computational domain, we find the shock front to propagate at slightly 
too large a pace. This artifact ceases with increasing numerical resolution, but 
future improvements of this point would be desirable. We have further deter- 
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mined the convergence rate of the method in numerical experiments and find 
it first-order accurate when shocks are involved and second-order accurate for 
smooth fiows. 
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